! Vind caculation program for RCR
! include "mod_param.f90"

! program vind4rcr
!     use mod_param 
!     implicit none 

!     real :: ctt, lamff 

!     call sol_ct4ff(ctt, lamff)
!     CT = T/(rho*vtip**2*Arot)
!     print *, ctt, lamff  
!     print *,CT

! end program vind4rcr

subroutine sol_ct4vf(thtot, ctt, lamh)
    use mod_param
    implicit none 

    real :: thtot
    real :: ctt, lamh ! output
    real :: a, &
            lam_n0, lam_n1, ct_n0, ct_n1 !tmp
    
    ! T = 5500

    ! Omega = 35.0
    ! R = 5.49
    ! chord = 0.365
    ! sigma = 0.0635
    ! Arot = pi*R**2
    ! cl_a = 6.0
    ! as = -3.0/180*pi 

    ! vtip = Omega*R 
    ! CT = T/(rho*vtip**2*Arot)

    ! eta0 = 50.0/180*pi 
    ! thtw = -10.0/180*pi
    ! th75 = eta0+0.75*thtw
    
    a = chord 
    ct_n0 = CT*2
    lam_n0 = sqrt(CT/2)

    do while(.true.)
        ct_n1 = 0.5*sigma*a*(thtot/3-lam_n0/2)
        lam_n1 = sqrt(ct_n1/2)
        ! print *, lam_n1
        if(abs(ct_n1-ct_n0)<=1e-6 .and. abs(lam_n1-lam_n0)<=1e-6) then 
            exit 
        else 
            ct_n0 = ct_n1 
            lam_n0 = lam_n1 
        end if 
    end do 

    ctt = ct_n1 
    lamh = lam_n1 
    print *, "VF: ", lamh*vtip

    return 
end subroutine sol_ct4vf

subroutine sol_ct4ff(thtot, ctt, lamff)
    use mod_param
    implicit none
    
    real :: thtot 
    real :: ctt, lamff ! output 
    real :: lam_n0, lam_n1, lam_in, &
            ctth, lamh, muz, ai, &
            frac  ! tmp

    ! T = 5500
    ! ! vair = 20

    ! Omega = 35.0
    ! R = 5.49
    ! chord = 0.365
    ! sigma = 0.0635
    ! Arot = pi*R**2
    ! cl_a = 6.0
    ! as = 3.0/180*pi 

    ! vtip = Omega*R 
    ! CT = T/(rho*vtip**2*Arot)
    ! mu = vair*cos(as)/vtip
    ! lam0 = vair*sin(as)/vtip

    ai = -as 
    muz = vair*sin(ai)/vtip

    frac = 0.5 ! fraction

    call sol_ct4vf(thtot, ctth, lamh)
    lam_n0 = lamh**2/sqrt((lamh+muz)**2+mu**2)+muz
    do while(.true.)
        lam_in = lamh**2/sqrt(lam_n0**2+mu**2)
        lam_n1 = lam_n0-(lam_n0-muz-lam_in)/(1+lam_in*lam_n0/(lam_n0**2+mu**2))*frac
        ! print *, lam_n1
        if(abs(lam_n1-lam_n0)<=1e-6) then 
            exit 
        else 
            lam_n0 = lam_n1 
        end if 
    end do 
    
    lamff = lam_n1 
    print *, "FF: ",lamff*vtip, "muz:", muz*vtip
    ctt = 2*sqrt(mu**2+lamff**2)*(lamff-muz)
    print *, "FF-CT: ", ctt 

    return 
end 

subroutine sol_ct4int(thtot0, thtot, ctint, lamint)
    use mod_param
    implicit none 

    real :: thtot0, thtot !input
    real :: out !output 
    real :: lamint, ctint, &
            lamint0, ctint0, &
            ai, muz !tmp 

    ai = -as 
    muz = vair*sin(ai)/vtip

    if(abs(vair-0.5)<=1) then 
        call sol_ct4vf(thtot0, ctint0, lamint0)
        call sol_ct4vf(thtot, ctint, lamint)
    else 
        call sol_ct4ff(thtot0, ctint0, lamint0)
        call sol_ct4ff(thtot, ctint, lamint)
    end if 

    out = 2*lamint*sqrt(mu**2+(muz+lamint0+lamint)**2)

    return 
end subroutine sol_ct4int